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Abstract 

We  present  a  semi-Lagrangian  method  for  integrating  the  three-dimensional  incompressible  Navier- 
Stokes  equations.  We  develop  stable  schemes  of  second-order  accuracy  in  time  and  spectral  accuracy 
in  space.  Specifically,  we  employ  a  spectral  element  (Jacobi)  expansion  in  one  direction  and  Fourier 
collocation  in  the  other  two  directions.  We  demonstrate  exponential  convergence  for  this  method,  and 
investigate  the  non-monotonic  behavior  of  the  temporal  error  for  an  exact  three-dimensional  solution. 
We  also  present  direct  numerical  simulations  of  a  turbulent  channel-flow,  and  demonstrate  the  stability 
of  this  approach  even  for  marginal  resolution  unlike  its  Eulerian  counterpart. 
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1  Introduction 


High-order  methods,  especially  spectral  methods,  are  particularly  effective  in  direct  numerical  simulations 
(DNS)  of  turbulent  ffows.  However,  most  Navier-Stokes  implementations  involve  semi-implicit  time  integra¬ 
tion  that  requires  unreasonable  small  time  step  sizes.  For  example,  for  a  flow  corresponding  to  Reynolds 
number  of  10,  000,  the  maximum  allowable  time  step  can  be  at  least  one  order  of  magnitude  smaller  than  the 
temporal  Kolmogorov  scale  [21].  It  can  be  projected  that  in  high  Reynolds  number  DNS  there  is  an  uneven 
distribution  of  resolution  in  space  and  time,  with  the  smallest  spatial  scale  approximately  matched  but  with 
the  temporal  scale  over-resolved  by  several  orders  of  magnitude.  This  inefficiency  of  currently  employed 
semi-implicit  schemes  for  DNS  of  inhomogeneous  turbulence  has  been  recognized  before,  and  attempts  have 
been  made  to  employ  fully  implicit  schemes  [3].  However,  this  requires  Newton  iterations  and  non-symmetric 
solvers  that  render  the  overall  approach  inefficient. 

Progress  can  be  made  by  employing  semi-Lagrangian  time-discretization,  which  could  increase  signifi¬ 
cantly  the  maximum  allowable  time  step  while  maintaining  the  efficiency  of  symmetric  solvers.  This  approach 
has  been  introduced  in  the  beginning  of  the  1980s  [17],  and  the  basic  idea  is  to  discretize  the  Lagrangian 
derivative  of  the  solution  in  time  instead  of  the  Eulerian  derivative.  The  work  here  is  an  extension  of  the 
semi-Lagrangian  scheme  proposed  in  [21]  but  formulated  in  the  context  of  simulating  turbulent  flows. 

The  semi-Lagrangian  method  depends  strongly  on  the  spatial  discretization.  Specifically,  its  accuracy 
is  particularly  sensitive  to  the  method  of  backward-integration  of  the  characteristic  equation  as  well  as 
the  interpolation  scheme  to  evaluate  the  solution  at  departure  points.  This  has  been  shown  by  Falcone  & 
Ferretti  [5]  who  conducted  a  rigorous  analysis  of  the  stability  and  convergence  properties  of  semi-Lagrangian 
schemes  for  advection-diffusion  equations.  It  has  also  been  shown  that  the  simplest  semi-Lagrangian  scheme 
with  linear  interpolation  is  equivalent  to  the  classical  first-order  upwinding  scheme  [15],  which  is  excessively 
dissipative  (see  [17]  and  [18]).  A  popular  and  effective  choice  for  interpolation  methods  in  previous  works  has 
been  the  cubic  spline  methods  [10];  see  also  [2].  The  idea  of  introducing  high  order  characteristic  methods 
has  first  been  presented  in  [4]  and  it  has  been  extended  into  the  spectral  frame  in  [9,  7]. 

An  intriguing  finding  is  that  the  error  of  semi-Lagrangian  schemes  in  solving  advection-diffusion  equations 
decreases  as  the  time  step  increases  in  a  certain  range  of  parameters,  and  this  has  initially  led  to  some 
erroneous  justifications  [13,  14].  The  analysis  in  [5]  showed  that  the  overall  error  of  the  semi-Lagrangian 
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method  is  indeed  not  monotonic  with  respect  to  time  step  At,  and,  in  particular,  it  has  the  form 


where  k  refers  to  the  order  of  backward  time  integration  and  P  to  the  (spatial)  interpolation  order;  similar 
conclusions  had  been  reached  earlier  in  [12].  Another  interesting  result  was  obtained  by  Giraldo  [6],  who 
combined  the  semi-Lagrangian  approach  with  a  spectral  element  discretization  for  the  advection-diffusion 
equation.  He  found  that  for  polynomial  order  P  >  4  the  combined  scheme  exhibits  neither  dissipation  nor 
dispersion  errors. 

The  extension  of  the  semi-Lagrangian  method  to  the  solution  of  Navier-Stokes  equations  was  presented 
in  the  pioneering  work  of  Pironneau  (1982)  [16].  He  demonstrated  nonlinear  stability  of  the  method  even 
as  the  viscosity  approaches  zero.  He  also  obtained  suboptimal  error  estimates,  which  were  improved  later 
by  Siili  (1988)  [19].  Most  of  the  previous  analysis  and  numerical  implementations  in  CFD  applications 
have  employed  the  Taylor-Hood  finite  element  and  are  first-order  in  time.  In  a  more  recent  paper  [1], 
an  error  analysis  was  conducted  for  the  fractional-step  method  for  incompressible  Navier-Stokes  equations. 
In  particular,  the  pressure-correction  version  of  the  fractional  scheme  with  first-order  time-stepping  was 
analyzed  and  an  extension  to  second-order  was  proposed  but  not  analyzed. 

In  this  paper,  we  present  a  semi-Lagrangian  method  for  simulating  three-dimensional  incompressible 
turbulent  flows  specifically  in  channel  domains,  extending  the  work  in  [21].  In  particular,  we  apply  a  Jacobi- 
based  spectral  element  discretization  along  the  inhomogeneous  direction  [8]  and  Fourier  collocation  along  the 
other  two  homogeneous  directions,  similar  to  [20] .  We  study  in  detail  the  dependence  of  the  overall  accuracy 
upon  the  time  step  for  an  exact  solution  of  the  three-dimensional  Navier-Stokes  equations.  We  then  present 
results  from  a  DNS  of  turbulence  for  48^  resolution  that  demonstrate  the  effectiveness  of  the  method. 

2  Formulation 

2.1  Advection-Diffusion  Equation 

Let  us  first  consider  an  advection-diffusion  equation  written  in  Eulerian  form 

=  (1) 
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and  in  semi-Lagrangian  form 


(2) 


Unlike  Lagrangian  formulations,  in  the  semi-Lagrangian  formulation  the  computational  mesh  is  fixed.  At 
each  time  step,  a  discrete  set  of  particles  arriving  at  the  grid  points  is  tracked  backwards  over  a  single  time 
step  along  its  characteristic  line  up  to  its  departure  points.  The  solution  values  at  the  departure  points 
are  then  obtained  by  interpolation.  For  example,  the  second-order  Crank-Nicolson  scheme  for  the  above 
equation  is 


At 


2A  =  u{x,  t),  =  Xa- 


dx 

dt 


(3) 

(4) 


Here  (j>2  denotes  the  value  of  (f)  at  the  departure  points  Xd  at  time  level  n,  and  Xa  is  the  position  of  the 
arrival  points  which  coincide  with  the  grid  points.  The  characteristic  equation  (4)  is  solved  backward,  i.e., 
we  solve  for  the  departure  point  at  time  level  n,  a;^,  with  the  initial  condition  =  Xa- 

The  departure  points  do  not  coincide  with  the  grid  points,  and  thus  a  search-interpolation  procedure  is 
needed.  Also,  the  overall  accuracy  and  efficiency  of  the  semi-Lagrangian  method  depends  critically  on  both 
the  accuracy  of  backward  integration  as  well  as  the  accuracy  of  the  interpolation  method.  In  the  following, 
we  provide  some  details  on  how  to  implement  both  algorithms. 

2.2  Backward  Integration 

We  solve  equation  (4)  for  one  single  time  step  in  order  to  obtain  Xd  =  a;(t")  by  the  explicit  second-order 
mid-point  rule 


At 

X  =  Xa - —U{Xa,t  ), 

.  „  At, 

Xd  =  Xa-Atu[X,t  +-^)- 


(5) 

(6) 


By  defining 


Ck,  —  Xd  X(^j 


we  can  re-write  the  explicit  mid-point  rule  as 


a  =  Atu{Xa  -  ^■u(Xa,t"),t"  +  ^). 


(7) 
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Similarly,  we  employ  implicit  integration  for  equation  (4)  setting 


x  =  Xa-  —  u{x,  +  — ), 

to  arrive  at  the  implicit  mid-point  rule 

.  ,  c(  „  At 

a  =  Atu(xa-  -,t  +-^)-  (8) 

This  is  the  backward-integration  algorithm  used  in  most  of  previous  semi-Lagrangian  schemes.  Although 
the  explicit  and  implicit  schemes  are  formally  of  second-order,  a  small  accuracy  improvement  has  been 
reported  for  the  implicit  scheme.  Equation  (8)  has  to  be  solved  iteratively,  but  numerical  experiments  show 
that  only  a  few  iterations  are  needed  for  convergence  (typically  around  five).  For  an  advection-diffusion 
equation  with  the  velocity  field  known  analytically,  the  additional  cost  associated  with  the  iterations  is 
negligible.  However,  for  a  velocity  field  known  only  in  numerical  form,  the  iteration  process  is  costly  because 
each  substep  requires  a  search-interpolation  procedure.  Our  numerical  results  show  that  the  two  methods  give 
almost  identical  results  in  practice,  and  for  more  general  problems,  especially  for  Navier-Stokes  equations, 
the  explicit  method  is  preferred. 

2.3  Search-Interpolation  Procedure 

We  consider  a  three-dimensional  channel  domain  with  spectral  elements//ip  along  the  inhomogeneous  di¬ 
rection  and  Fourier  collocation  along  the  two  homogeneous  directions  (streamwise  and  spanwise).  In  its 
current  implementation,  we  first  locate  in  which  spectral  element  the  departure  point  resides  and  perform 
interpolation  in  modal  space  along  all  three-directions.  This  involves  the  entire  Fourier  expansions  which  is 
computationally  expensive.  Currently,  we  are  working  on  an  implementation  that  employs  a  more  localized 
interpolation  in  the  Fourier  directions  to  reduce  that  cost. 

2.4  Incompressible  Navier-Stokes  Equations 

We  consider  the  incompressible  Navier-Stokes  equations  in  Lagrangian  form 

du  o  ,  ■. 

—  = —\7p  +  V  ■  u  =  0.  (9) 

dt 
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We  employ  a  stiffly-stable  scheme  to  discretize  the  above  equations  [8].  A  second-order  time-discretization 


corresponds  to 


=  (_Vp  -H 


where  is  the  velocity  u  at  the  departure  point  at  time  level  and  ^  is  the  velocity  at  the  departure 
point  x'^~^  at  time  level  The  departure  point  is  obtained  by  solving 

at 


and  also 


The  point  a;^  ^  is  obtained  by  solving 


^+i/2  =  3/2u"-l/2r 


=  u"(x,t),  a;(t"+^)  =  Xa- 


By  using  the  above  characteristic  equations,  the  resulting  scheme  is  second-order  accurate  in  time. 

Specifically,  for  computational  convenience  we  use  the  following  three  substeps  to  solve  equation  (10) 

u-2u]t  + 

a  2  d  _  n  /'I 


|u"+i  -  ft 


=  -Vp" 


The  discrete  divergence-free  condition  results  in  a  consistent  Poisson  equation  for  the  pressure,  i.e. 


V =  ttV  •  *, 


with  accurate  pressure  boundary  conditions  of  the  form  [8] 


=  -i^-[u  +  W  X  u;"+y 


where  n  is  the  unit  normal,  and  lu  is  the  vorticity. 
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3  Numerical  Results 

3.1  Convergence  Rate 

We  first  present  results  from  comparisons  with  an  exact  three-dimensional  solution  to  the  incompressible 
Navier-Stokes  equations,  given  by 

u  =  sin(mx)  cos{ly)  cos{nz)e~*^^^ 

m  +  n  .  \  ■  n  \  /  \  -tiRe 

V  = - j — cos{mx)sin{ly)cos{nz)e  ' 

w  =  cos{mx)  cos{ly)  sin(nz)e“‘/^® 

where  Re  is  the  Reynolds  number  and  m,l,n  define  the  wavenumbers  along  the  three  directions.  We 
determine  the  pressure  p{x,  y,  z,  t)  from  the  Navier-Stokes  equations  (assuming  that  no  forcing  is  applied  in 
the  y-direction,  i.e.  fy  =  0)  to  be 

p{x,  y,  z,  t)  =  —  ^  (w^  +  1"^  +  11?  —  1)  cosimx)  cos( ly)  cos(nz)e~*''^‘^  + 

VRe 

TfiiTTl  ~\~  77-)  .9/  \  /c-\i  \  2/  \  _ 2// Rp 

- — - ^  Sin^(777,x)  COs(2/y)  COS^(77-2)e  + 

4/^ 

cos^(mx)  cos(2/y)  cos^(nz)e“^*/'^'^  -I- 
cos^{mx)  cos{2ly)  sii?{nz)e~^*^^^ . 

With  this  expression  for  pressure  we  can  now  evaluate  the  extra  forces  along  the  two  other  directions  fx 
and  fz,  which  are  computed  so  that  the  above  is  an  exact  Navier-Stokes  solution.  The  domain  used  was  a 
cube  of  size  27r^. 

We  first  tested  that  exponential  convergence  is  realized  with  errors  almost  identical  to  the  Eulerian 
approach.  In  figure  1  we  plot  the  L2  error  for  all  three  components  of  velocity  for  wavenumbers  m  =  I  =  n  = 
1.  Only  one  element  was  used  along  the  y-direction  for  this  solution,  with  P  denoting  the  Jacobi  polynomial 
order.  Similar  results  hold  for  a  multi-element  discretization.  Also,  the  number  of  collocation  points  along 
the  two  Fourier  directions  is  set  to  M  =  N  =  P  for  this  case.  The  Reynolds  number  was  set  to  Re  =  1  and 
the  final  time  of  integration  was  T  =  1 .  The  non-monotonic  behavior  of  the  temporal  error  obtained  by  the 
theoretical  analysis  for  advection-diffusion  equations  is 

A 

OiA?  +  ^), 
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Figure  1:  Errors  versus  polynomial  order  for  the  exact  three-dimensional  solution  {Re  =  1). 

where  As  is  an  equivalent  grid  spacing.  This  behavior  is  realized  in  our  computations  of  the  exact  Navier- 
Stokes  solutions,  as  shown  in  figure  2.  We  note  that  the  u  and  w  velocity  components  correspond  to  identical 
curves  but  the  normal  velocity  component  v  exhibits  a  different  behavior.  This  different  behavior  can  be 
attributed  to  the  fractional  stepping  scheme  and  it  is  similar  to  the  Eulerian  approach.  If  we  increase  the 
resolution  along  the  y-direction  only  from  P  =  6  to  P  =  8  (figure  3),  we  observe  a  large  reduction  in  the 
error  for  all  components  but  also  a  different  qualitative  trend.  This  is  consistent  with  the  aforementioned 
error  estimate  and  the  fact  that  exponential  convergence  is  achieved.  The  results  in  both  plots  suggest  that 
both  error  terms  are  comparable  for  this  resolution. 

At  higher  resolution,  the  first  term  in  the  error  dominates  and  this  behavior  is  similar  to  the  one  obtained 
in  the  Eulerian  approach.  A  comparison  of  the  two  approaches  is  shown  in  figure  4;  for  the  semi-Lagrangian 
scheme  the  P2  error  is  slightly  larger  compared  to  the  Eulerian  scheme.  We  also  see  from  figure  4  that  a 
second-order  time-accuracy  is  obtained.  In  figure  5  we  plot  the  errors  for  Reynolds  number  Re  =  10®  and 
final  time  of  integration  T  =  1.  We  see  that  the  same  non- monotonic  behavior  is  obtained  as  in  the  low 
Reynolds  number  case  Re  =  1  studied  in  the  previous  figures. 
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Figure  2:  Errors  versus  time  step  for  the  exact  three-dimensional  solution;  resolution  P  =  6  =  M  =  N. 
{Re=  1). 
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Figure  3:  Errors  versus  time  step  for  the  exact  three-dimensional  solution;  resolution  P  =  8  and  M  =  N  =  6. 
{Re  =  1). 
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Figure  4:  Comparison  of  Eulerian  (EU)  and  semi-Lagrangian  schemes  (SL).  Errors  versus  time  step  for  the 
exact  three-dimensional  solution;  resolution  P  =  M  =  N  =  12.  {Re  =  1). 


Figure  5:  Reynolds  number  Re  =  10®.  Errors  versus  time  step  for  the  exact  three-dimensional  solution; 
resolution  P  =  M  =  N  =  12. 
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3.2  DNS  of  Turbulent  Channel-Flow 


Here  we  perform  direct  numerical  simulations  of  turbulent  channel-flow  in  a  periodic  domain  (in  x-  and 
z-directions)  of  size  (x,y,z):  27r  x  2  x  27r.  The  Reynolds  number  based  on  the  half-channel  width  and  the  wall 
shear  velocity  is  i?e*  =  116.  We  interpolated  an  initial  turbulent  64^  held  to  a  48^  held  and  ran  the  simulation 
for  several  eddy  turn-over  times  (50  convective  units)  with  the  de-aliased  Eulerian  code.  We  note  here  that 
the  corresponding  Eulerian  simulation  without  de-aliasing  was  unstable  even  at  very  small  time  steps.  We 
then  collected  statistics  for  T  =  5  and  T  =  20  convective  units  for  the  Eulerian  simulation;  a  time  step  of 
At  =  1/200  was  employed.  The  semi-Lagrangian  simulation  was  performed  with  time  step  At  =  1/20,  and 
statistics  were  gathered  for  total  time  T  =  20  convective  units.  Of  course,  no  de-aliasing  is  required  in  this 
case.  In  wall  units  the  time  step  used  in  the  semi-Lagrangian  method  is  At+  =  At  ■  u^/i^  =  0.2645,  which  is 
smaller  than  the  Kolmogorov  time  scale;  Ur  is  the  wall  shear  velocity.  The  Kolmogorov  scale  was  estimated 
at  =  5  to  be  «  0.33  using  an  approximate  value  for  the  dissipation  rate  e  «  0.10  [11].  We 

note  here  that  the  time  Kolmogorov  scale  in  [3]  was  overestimated;  the  value  of  2.4  (in  wall  units)  reported 
would  correspond  to  a  value  of  dissipation  rate  of  e  «  0.002  compared  to  about  0.16  at  the  wall  [11]. 

A  comparison  of  the  mean  velocity  profile  is  shown  in  figure  6;  no  significant  differences  are  observed 
despite  the  very  large  time  step  employed  in  the  semi-Lagrangian  simulation.  However,  some  differences  are 
present  in  figure  7,  where  we  plot  the  rms  values  of  all  three  velocity  components  for  both  approaches.  The 
more  pronounced  differences  correspond  to  the  streamwise  turbulent  intensity,  which  mat  be  associated  with 
time-averaging  errors.  To  evaluate  this  we  compare  these  profiles  obtained  by  averaging  also  over  T  =  5  and 
we,  indeed,  obtain  noticeable  differences  as  shown  in  figure  8. 

4  Summary 

We  have  developed  a  spectral  semi-Lagrangian  algorithm  for  simulating  turbulent  channel  flow,  and  have 
demonstrated  its  accuracy  and  its  stability.  The  overall  error  is  comparable  to  the  Eulerian  scheme  but 
the  stability  is  greatly  enhanced.  However,  in  its  current  implementation  the  method  employs  global  in¬ 
terpolation,  which  makes  it  computationally  very  expensive.  Specifically,  the  method  is  approximately  ten 
times  slower  than  the  Eulerian  method  for  the  48^  simulation  presented,  so  with  the  gain  factor  of  ten  (for 
the  semi-Lagrangian)  in  the  time  step,  the  overall  simulation  cost  is  the  same  for  both  methods.  Better 
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Figure  6:  Profile  of  mean  velocity  from  the  48^  DNS.  The  solid  line  corresponds  to  the  Eulerian  simulation 
and  the  dash  line  to  the  semi-Lagrangian  simulation.  Reynolds  number  i?*  =  116. 


Figure  7:  Comparison  of  turbulence  intensities  from  the  48^  DNS.  The  solid  line  corresponds  to  the  Eulerian 
simulation  and  the  dash  line  to  the  semi-Lagrangian  simulation.  Reynolds  number  i?*  =  116. 
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Figure  8:  Comparison  of  turbulence  intensities  from  the  Eulerian  48^  DNS.  The  solid  line  corresponds  to 
averaging  over  T  =  40  and  the  dash  line  to  averaging  over  T  =  5;  Reynolds  number  i?*  =  116. 


local  interpolation  procedures  need  to  be  implemented  that  do  not  compromise  accuracy  while  providing  a 
speed-up  factor  that  will  make  this  method  more  efficient  than  the  Eulerian  approach  for  DNS  of  turbulence. 
This  has  been  done  in  the  context  of  fully  three-dimensional  geometries,  where  hexahedra  and  tetrahedra 
spectral/ft-p  elements  form  a  natural  framework  for  local  Lagrangian  interpolation  [21].  In  that  case  a  speed¬ 
up  factor  of  about  five  (on  average)  was  achieved  in  favor  of  the  semi-Lagrangian  method.  A  similar  domain 
decomposition  scheme  should  also  be  adopted  for  the  channel  geometry;  we  will  report  such  work  as  well  as 
practical  high  Reynolds  number  DNS  in  future  publications. 
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